• Review

• Multiple factors

• Interactions

The linear model estimates the mean of the outcome given the covariates.

• the mean body weight for a given heart girth.

• the mean body weight given heart girth, umbilical girth, length and sex.

• the mean right wing length for different areas.

It also gives uncertainties, i.e. confidence intervals.

• For factor (grouping) variables, it estimates a separate mean for each level.

• The first level (baselines) goes into the intercept.

• Each other level is represented by the difference to the baseline.

• The overall F-test tells us whether the grouping variable is important.

\[ F = \frac{\mathsf{Variation\ explained}}{\mathsf{Variation\ unexplained}} = \frac{(\sigma^2_\mathsf{total} - \sigma^2_\mathsf{res})/p}{\sigma^2_\mathsf{res}/(n-p-1)} \]

• The numerator is variation explained by the model. In this case this is variation between groups, as the model is just giving a different mean to each group.

• The denominator is residual variation, which in this case is variation within groups (residuals are values minus fit, which is the group mean).

• Thus F is a ratio of between-group variation to within-group variation.

• The overall F-test’s P-value suggests average right wing length differs between areas.

• Summary table suggests that Area 5 has largest birds, having 4.14mm larger ring wing lengths compared to those in Area 1.

• Birds in Area 6 are smallest, having 9.27 mm smaller right wing lengths compared to those in Area 1.

• We can get confidence intervals for the mean using the predict function

new_data <- data.frame(Area=factor(1:6))
predict(mod, new_data, interval="confidence")

       fit      lwr      upr
1 389.6000 388.5335 390.6665
2 383.5964 382.4229 384.7699
3 384.3750 379.9940 388.7560
4 382.9107 380.5690 385.2525
5 393.7462 392.2092 395.2831
6 380.3333 375.2746 385.3921

library(visreg)
visreg(mod)

• We include a block of indicator variables for each factor.

• The first level, or baseline of each of the factors is what contributes to the intercept.

• Subsequent levels are treated as differences compared to the baseline level.

• The overall F-test and P-value for the model is testing whether any of the factors in the model are important.

• We need another way to compute P-values for each factor separately.

• Fit linear model.

• Do anova() to check the factor is important (small P).

• If it is important, look at the summary() to see how the levels differ (or visreg() to visualise them).

• Multiple testing problem on the latter, which is why we do the overall test first.

• Order can be important in the anova table.

• Each row’s P-values are adjusting for previous rows.

anova(lm(R.Wing.Lth ~ Sex + Area, data=petrels))

Analysis of Variance Table

Response: R.Wing.Lth
           Df Sum Sq Mean Sq F value    Pr(>F)    
Sex         1   1107 1106.58  14.328 0.0001668 ***
Area        5  11224 2244.84  29.066 < 2.2e-16 ***
Residuals 696  53754   77.23                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(lm(R.Wing.Lth ~ Area + Sex, data=petrels))

Analysis of Variance Table

Response: R.Wing.Lth
           Df Sum Sq Mean Sq F value  Pr(>F)    
Area        5  11870 2374.10 30.7393 < 2e-16 ***
Sex         1    460  460.29  5.9597 0.01488 *  
Residuals 696  53754   77.23                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

• Each row represents a single variable (factor or numeric).

• Each row adjusted for what is above it (but not below).

• Order can be important but probably tells you something if it changes things.

• If a factor is important, look at summary table to see which groups differ (or visualise model).

• From ANOVA table, both Area and Sex are important to the right wing length.

• From summary table, we see that males have 2.04mm larger right wing length than females after accounting for differences between areas.

• This doesn’t mean that male birds have a 2.04mm larger right wing length than females.

• In fact

  tapply(petrels$R.Wing.Lth, petrels$Sex, mean, na.rm=TRUE)

    Female     Male 
    385.2222 388.3309 
    

  they’re 3.11mm larger across all areas.

• Similarly, birds in area 5 have 3.76mm larger right wing length than those in area 1 after accounting for differences in sex.

• Again,

  tapply(petrels$R.Wing.Lth, petrels$Area, mean, na.rm=TRUE)

         1        2        3        4        5        6 
    389.6000 383.5964 384.3750 382.9107 393.7462 380.3333 
    

  they’re 4.15mm larger across both sexes.

• Always look at the anova table first, particularly if you have factors.

• Look at the summary table next, and interpret the estimated coefficients there for the significant variables from the anova table.

• Use visreg to visualise the effects - easier to see effect sizes/uncertainty.

• Remember that ‘not significant’ does not mean ‘not important’.

• Our current models have been assuming that each effect is additive.

• e.g. the difference in the right wing length between the sexes is common across all areas, and similarly the area effects are common for males and females.

• This is often not the case.

• Suppose we have some medical condition and we have a numerical measure describing how good/bad a patient is.

• Suppose further that we have two drugs A and B to treat that condition, and that we have the option of combining them.

• Then there are 4 possible treatments: No drugs, just drug A, just drug B, or both drugs.

• Our current linear models, would assume that we can add the effect of drug A and B together for the “both drugs” case.

The model equation would be \[ \mathsf{mean}(y) = \alpha + \beta_A z_A + \beta_B z_B \] where \(z_A\) and \(z_B\) are indicators for whether drug A or B are included, giving

• For the combined drugs treatment, we have \(\mathsf{mean}(y)=\alpha + \beta_A + \beta_B\)

• We’re thus assuming that the individual benefit of each drug add together when combined.

• But, it might be that the combined effect of the drugs gives is smaller or larger than the sum of individual effects.

• We need a bit more flexibility in the model.

We add another indicator variable \(z_{AB}\) which is 1 when both drugs are in use. \[ \mathsf{mean}(y) = \alpha + \beta_A z_A + \beta_B z_B + \beta_{AB} z_{AB} \]

Mathematical aside: the new variable is the product of the existing ones, as \(z_{AB}=1\) only if both \(z_A=1\) and \(z_B=1\). \(z_{A} \times z_{B}\) has this property.

The means for the neither drug or single drug treatments are still the same, as \(z_{AB}=0\) in those cases. The mean in the both drugs treatment is the only one that changes to \[ \mathsf{mean}(y) = \alpha + \beta_A + \beta_B + \beta_{AB}. \] Thus, \(\beta_{AB}\) can be used to measure whether or not the effects of the two drugs are additive or not.

• \(\beta_{AB}\) can be used to measure whether or not the effects of the two drugs are additive or not.

• An alternate way to think about it is that the effect of drug B may differ depending on whether the patient is already given drug A.

• e.g. the effect for drug B may be to increase the diagnostic measurement by 5 units in the absence of other treatment, but if the patient is on drug A as well, then the effect of drug B may be less pronounced.

mod3 = lm(R.Wing.Lth ~ Sex + Area + Sex:Area, data=petrels)
anova(mod3)

Analysis of Variance Table

Response: R.Wing.Lth
           Df Sum Sq Mean Sq F value    Pr(>F)    
Sex         1   1107 1106.58 14.3035  0.000169 ***
Area        5  11224 2244.84 29.0165 < 2.2e-16 ***
Sex:Area    5    296   59.15  0.7646  0.575486    
Residuals 691  53459   77.36                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

As the interaction here is between factors, it should be evaluated using the ANOVA table.

• From the anova table we see the interaction term is not significant (P=0.575).

• Thus, there is no evidence in the data to suggest that the comparative size of male and female petrels differ between areas.

• While the interaction term isn’t significant (thus all the interaction coefficients could be zero in the population) let’s assess what each row represents.

• The Intercept contains the first levels of each factor, so females in area 1 have on average 389.12mm right wing lengths.

• The SexMale term is the difference to the baseline, so this describes how males in area 1 differ. They have right wing lengths 0.8mm larger.

• The Area2 term is the difference to the baseline, so this describes how females in area 2 differ from those in area 1. Their right wing lengths are 8.37mm smaller.

• The SexMale:Area2 is the differential effect of males in area 2 over and above the SexMale and Area2 effects. Compared with females in area 1, males in area 2 are \(0.8+-8.37+2.73=-4.84\)mm larger.

• When you have interactions, visualise the effect of one variable within the levels of another.

• Allows you to see how effects interact.

• Better way to present results from the model than summary table.

library(visreg)
visreg(mod3, "Area", by="Sex")
visreg(mod3, "Area", by="Sex", overlay=TRUE)

Consider \[ \mathsf{mean}(y) = \alpha + \beta_1 x + \beta_2 z \] where \(x\) is a numeric variable and \(z\) is an indicator variable for the 2nd level of a factor variable.

Then if we’re in the first level of the factor variable, \(z=0\) so the equation is \[ \mathsf{mean}(y) = \alpha + \beta_1 x \] If we’re in the second level, \(z=1\) so that \[ \mathsf{mean}(y) = \alpha + \beta_1 x + \beta_2 = (\alpha + \beta_2) + \beta_1 x. \]

The model fit is thus two parallel lines, separated by \(\beta_2\).

c1 = lm(Weight ~ BirthWeight + Breed + Treatment + Age, data=calf)
anova(c1)

Analysis of Variance Table

Response: Weight
             Df Sum Sq Mean Sq  F value    Pr(>F)    
BirthWeight   1  10952   10952  390.912 < 2.2e-16 ***
Breed         2    888     444   15.853 2.037e-07 ***
Treatment     1    657     657   23.438 1.685e-06 ***
Age           1 212268  212268 7576.597 < 2.2e-16 ***
Residuals   542  15185      28                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

visreg(c1, "Age", by="Treatment", overlay=TRUE, gg=TRUE)

visreg(c1, "Age", by="Breed", overlay=TRUE, gg=TRUE)

We could add another term using the product of the numeric variable \(x\) and the indicator for the second level of the factor \(z\) \[ \mathsf{mean}(y) = \alpha + \beta_1 x + \beta_2 z + \beta_3 x z. \]

Then if we’re in the first level of the factor variable, \(z=0\) so the equation is \[ \mathsf{mean}(y) = \alpha + \beta_1 x \] whereas if we’re in the second level, \(z=1\) so that \[ \mathsf{mean}(y) = \alpha + \beta_1 x + \beta_2 + \beta_3 x = (\alpha + \beta_2) + (\beta_1 + \beta_3) x. \] The model fit is thus two separate lines.

c2 <- lm(Weight ~ BirthWeight + Breed + Treatment + Age +
           Age:Treatment + Age:Breed, data=calf)
anova(c2)

Analysis of Variance Table

Response: Weight
               Df Sum Sq Mean Sq  F value    Pr(>F)    
BirthWeight     1  10952   10952  420.039 < 2.2e-16 ***
Breed           2    888     444   17.034 6.707e-08 ***
Treatment       1    657     657   25.184 7.088e-07 ***
Age             1 212268  212268 8141.130 < 2.2e-16 ***
Treatment:Age   1    494     494   18.936 1.617e-05 ***
Breed:Age       2    637     319   12.224 6.427e-06 ***
Residuals     539  14054      26                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

c3 <- lm(log(Weight) ~ log(BirthWeight) + Breed + Treatment + Age + 
         Age:Treatment + Age:Breed, data=calf)
anova(c3)

Analysis of Variance Table

Response: log(Weight)
                  Df Sum Sq Mean Sq   F value    Pr(>F)    
log(BirthWeight)   1  3.613   3.613  599.9067 < 2.2e-16 ***
Breed              2  0.182   0.091   15.0937 4.187e-07 ***
Treatment          1  0.113   0.113   18.8347 1.703e-05 ***
Age                1 58.253  58.253 9672.8215 < 2.2e-16 ***
Treatment:Age      1  0.086   0.086   14.2626 0.0001767 ***
Breed:Age          2  0.070   0.035    5.8391 0.0030987 ** 
Residuals        539  3.246   0.006                        
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

• We know here that we have multiple calves measured through time.

• Residuals from the same calf likely more similar to each other than to those of other calves.

• The residuals are unlikely to be independent.

• In this case can check by plotting the residuals versus the calf identifier.

calf_pred = broom::augment(c3, data=calf)
ggplot(calf_pred, aes(x=Calf, group=Calf, y=.resid)) + geom_boxplot() + ylab("Residual")

• Review

• Multiple factors

• Interactions